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0^ ' We present a method of simulating lattice QCD at nonzero chemical po- 



tential in the chiral limit. By adding a weak four-fermi interaction to the 
standard staggered fermion SU(3) QCD action, we produce an algorithm in 
which the limit of massless fermions is well-behaved and physical. Using con- 



^ ' figurations at zero chemical potential, and an exact fugacity expansion of the 



fermion determinant, we can simulate QCD at nonzero chemical potential and 
evade the notorious problem of the complex action. Small lattice simulations 
give physical results : At strong gauge coupling the critical chemical poten- 
tial agrees with theoretical expectations and at weak gauge coupling is 
nonzero in the low temperature confined phase of QCD and jumps to zero in 
the high temperature quark-gluon plasma phase. In all these simulations the 

quarks are exactly massless and there is a Goldstone pion. 
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I. INTRODUCTION 



One of the goals of Lattice Gauge Theory is a thorough investigation of the phase diagram 
of QCD in the temperature-chemical potential plane. Considerable progress has been made 
in simulating lattice QCD at nonzero temperature but there has been little progress 
in understanding the theory at nonzero chemical potential n even though a sound lattice 
formulation was presented twelve years ago 0. One expects on the basis of elementary 
physics arguments that as the chemical potential /i increases QCD will experience a 
transition to a quark-gluon plasma at one-third the mass of the proton, and that the chiral 
condensate will serve as an order parameter. The basic difficulty in simulating QCD at finite 
density and investigating the transition quantitatively is that the effective action resulting 
from the Grassmann integration over the fermions is complex due to the introduction of 
the chemical potential in the Dirac matrix. The standard simulation algorithms 03 
lattice QCD with dynamical fermions do not apply to this situation. In order to avoid 
the determinant, early simulations used the quenched approximation, but they led to an 
unphysical value of the critical chemical potential /ic, namely, its proportionality to the pion 
mass [^]. The physical and mathematical reasons for this failure have been the subject of 
considerable debate [^]. Recent work by M. Stephanov [|^, however, has shown that the 
quenched model is the rij — > limit of a theory with both quarks and conjugate quarks, and 
is not relevant to QCD with dynamical fermions. 

Here we wish to advocate a new approach to simulating QCD with dynamical fermions 
at nonzero chemical potential which has two unique features. First, we approach the chiral 
limit of lattice QCD by adding a four-fermi term to its action but setting the bare masses 
of the quarks to zero. The four-fermi term is a perturbatively irrelevant interaction so it 
does not have an effect on the values of physical quantities in the theory's continuum limit. 
However, in the new algorithm, called xQCD the dynamical mass of the quark enters its 
propagator directly and makes the limit of massless quarks more accessible, as we will see 
below. The physics of massless Goldstone pions can be studied in xQCD without the need 



2 



for extrapolation. In addition, the algorithm runs one to two orders of magnitude faster 
than the traditional one because a physical quantity, the dynamical quark mass rather than 
the bare lattice mass, controls the convergence of the inversion of the quark propagator, the 
most time consuming step in the procedure. xQCD approaches the chiral limit of QCD in 
a different direction than the standard action and appears to avoid some of its problems. 
Second, the fermion determinant is calculated exactly at any /i and any bare quark mass 
rriq, given a configuration produced by the xQCD algorithm at vanishing /i and any other 
ruq. This can be done because of the simple way and fi enter the theory's action. A 
large ensemble of configurations will generally be needed to produce observables at nonzero 
fi values with acceptable statistical errors. Preliminary simulations of the method have been 
successful and will be presented below. 

II. FORMULATION 

Consider QCD with a chiral 4-fermion interaction, 'xQCD' 0]. The molecular dynamics 
Lagrangian for the lattice theory is 

L = - \Re{TTnUUUU)]+Y,i^^M^Mij 

□ -J s 

+^ E(^? + ^8 + m + ¥2 + w (1) 

where 

M=^ + mg + 4E(^^ + ^e7r,) (2) 

and P = 7 = l/G^, e = ( — 1)^+?'+^+*^ the 9i parametrize the SU(3) link variables U, s 
labels dual sites ^ on the lattice of volume V = rit x n^, and we have introduced auxiliary 
fields a and vr to formally linearize the 4-fermi term. We use staggered quarks in Eq. 1 
because of their relatively good chiral symmetry properties. For simplicity we consider a 
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theory where the 4-fermion operator has the 

f/(l) X U{1) C SU{Nf) X SU{Nf) 



flavor symmetry generated by (1,^75^5) |T^. This Lagrangian describes 8 flavors. For Nf 
which is not a multiple of 8 we use "noisy" fermions and multiply the fermion kinetic 
term by Nf/8. 

The Dirac operator of standard lattice QCD becomes singular as ruq 0. However, 
inspecting Eq. 2, we see that in xQCD the Dirac operator remains non-singular at rriq = 
because the auxiliary field a develops a vacuum expectation value due to chiral symmetry 
breaking. Conjugate gradient inversion of the regular QCD Dirac operator requires a number 
of iterations which diverges as V — > 00 and rrig — > 0. Inversion of the xQCD Dirac operator 
requires a finite number of iterations even at rUg = 0. In addition, the scale of the 'time' 
step in the molecular dynamics algorithm of xQCD is set by the dynamical quark mass and 
can be chosen several times larger here than in the ordinary lattice QCD algorithm, for the 
same systematic error or acceptance rate [Q. 

In order to circumvent the difficulty of investigating the finite-density transition for 



dynamical fermions, we build on the method of reference |TT|]. The method was inspired 



by the classic work of Yang and Lee ||12| who showed that the distribution of the zeros 
of a partition function determines the equation of state for a many-body system. For an 
Ising model in a magnetic field and for a lattice gas, it was shown that the zeros of the 
corresponding partition function lie on a unit circle in the complex fugacity plane. For 
a finite system these zeros will never lie in the physical range for the fugacity, namely the 
positive real axis. However, for a temperature below the phase transition, they will approach 
the real axis as the volume of the system tends to infinity. The singularities of the free energy 
in the thermodynamical limit are therefore obtained as the infinite-volume extrapolation of 
the zeros of the partition function. A detailed investigation of the finite-size scaling of the 
zeros enables one to evaluate the order of the phase transition and the critical exponent for 
the order parameter |]13|. 



This proposed study of lattice QCD at finite density is based on expanding the Grand 
Canonical Partition Function (GCPF) as a polynomial in the fugacity variable [e^"'^) where 
rir is the temporal extent of the lattice. 

The GCPF is given as an ensemble average of the determinant of the Dirac operator 
normalised with respect to the determinant at /i = 0: 

J[dU][dW] detM(/i) detM(^)e-^«[^'^^l 



J[dU][dW] det M(/i = 0) detM(/i = 0)e-^«[^'^^l 
The Dirac fermion matrices M and M are given by: 



(3) 



-2zM,.,(^) = F4 + G^y + V^ye^ + V^e'^ (4) 

where Y = 2i{mq + ^ I]<a;,i>('^(5^) + ^e7r(a;)))(5^j^, G contains all the spacelike links and V 

all the forward timelike links. 

The determinant of the Dirac operator is complex (for /i 7^ 0) in the Nc = ?> case because 
favors forward propagation through Eq. 4 while hinders backward propagation |0. 
The determinants of the fermion matrices M and M are related to that of the propagator 

matrix M 



P 



( -GV -YV 
y -V Oj 



(5) 



by det(2iM) = e^^'^'"'"* det(P-e-'^) and det(2zM) = e^=^"'"' det{{P-^y -e-^") on an nlxrit 
lattice. 

Since the eigenvalues of P have a Z{nt) symmetry (since P is proportional to V) and 
because the complex conjugate configuration is equally probable in the ensemble average, 
the GCPF can be expanded as 

Z= < hn\ > e"''"* = E e-(^"-"'^)/^. (6) 

n=-2Nc.ns^ n=-2Ncns^ 



The bn are determined from the eigenvalues of P 



nt 



For Nc = 3, the tunneling between the different Z{3) vacua should eliminate the triality 
non-zero coefficients. In the simulations described below we do observe strong cancellations 
in these coefficients for n = around zero. 

The zeros of this polynomial are the Lee- Yang zeros in the complex fugacity plane. Their 
distribution and behaviour with respect to finite volume scaling will be discussed in a future 
publication. 



III. RESULTS 

We consider a description of the system in terms of the canonical partition functions 
for fixed particle number [|15|. The chemical potential as a function of the baryon number 



density is obtained from the local derivative of the energy, e„, of the state with n fermions 
with respect to n. 

1 de 

where p is the fermion density, ^^^"^3 . 

Although some of the e„ are determined with large error and have imaginary part m 
because the corresponding 6„ is negative (in principle, many more measurements should 
give all triality zero 6„ positive), it is clear from the simulations that, for all mod(n,3)=0, 
their real parts form a continuous curve in n, as shown in Fig. 1. 

With this assumption we made a cubic spline fit to a randomly selected subset of 
of the 2n^ e„'s with triality zero and evaluated the derivative. This process was performed 
n^/3 times and an estimate of the fitting error and the mean estimated from the distribution 
of the corresponding /i(p) resulted. 
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FIG. 2. Fermion density vs. chemical potential, 4^ lattice, P = 0.5 
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FIG. 3. Fermion density vs. chemical potential, 6^ lattice, /3 = 5.0 

Our preliminary simulation results from small lattices are very encouraging. In Fig. 2 
we show strong gauge coupling (/5 = 0.5, = 0) results of the induced fermion number vs. 
chemical potential // on a 4^* lattice at various four-fermi couplings 7. Since the pion is strictly 
massless in these simulations and since [ic is distinctly nonzero for each parameter choice, 
the simulation is not afflicted by the diseases of the quenched standard action. In fact, for 
large 7 where the action reduces to the usual lattice QCD action, mean field analyses [jl6 



which should be reliable at strong gauge coupling, predict /ic near 0.60, in rough agreement 
with our figure. 

Next we simulated a weaker gauge coupling in order to approach the continuum limit of 
the theory and test the simulation method more strenuously. We chose (5 = 5.0, rrig = for 
various 7 on a larger 6^ lattice. Simulations at zero /i were performed and measurements 
of the chiral condensate and heavy quark potential indicated a finite 'temperature' quark 
gluon plasma transition between 7 of 5.0 and 6.0. (At 7 = 5.0 we measured < a >= .20 and 
the Wilson Line < WL >= .01, while at 7 = 6.0 we found < a >= .03 and < WL >= .28.) 
The task for our new algorithm was to find this transition and measure the critical chemical 
potentials in the 'cold' hadronic phase for 7 below 5.0 and do the same in the 'hot' plasma 
phase for 7 above 6.0. Our simulation results, shown in Fig. 3, show this transition very 
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clearly - for 7 of 5.0 and smaller the /ic is above 0.50, while at 7 = 6.0 /ic falls abruptly to 
zero, indicating the presence of unconfined, 'massless' quarks. Clearly, this estimate of /ic 
is crude. Simulations of simpler models on much larger lattices suggests that the 'tails' 
in the curves in Fig. 3 are finite size effects. The finite temperature transition of xQCD 
is being studied on larger lattices for 7 > 10 at variable (3 and more observables are being 
measured |^ to obtain quantitative predictions for pure QCD. 

Preliminary measurements, on 4^ and 6^ lattices, of the fermion energy density and 
condensate using the stochastic estimator method and including the ratio of the determinants 
of the Dirac operator at yU to that at /i = are consistent with the behaviour of the number 
density found by the above method. 

At the present time extension to larger lattices is limited by the algorithm used to 
determine the eigenvalues of the propagator matrix. The cpu time and storage, to a first 
approximation, scale as but are independent of rit. The algorithm has been successfully 
implemented on an 8^ lattice and should be viable on a 10^ x rit lattice at its present stage 
of development. 

IV. CONCLUSIONS 

Clearly much more work along these lines is called for - much larger lattices must be used 
to extract continuum physics. Although in principle the four fermi term is irrelevant and 
the value of 7 should not effect the values of observables in physical units in the continuum 
limit, practical simulations on finite lattices should be done with the four fermi coupling 
sufficiently small. This issue is discussed briefly in [§, but more quantitative work is needed 
and is underway. Nonetheless, the method has passed several nontrivial tests which we find 
encouraging. 

We are also hopeful that xQCD can be used to improve the action of lattice QCD, speed 
up spectroscopy and matrix element calculations of conventional hadronic phenomenology, 
and determine the universality class of the finite temperature hadronic matter/quark gluon 
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plasma transition 
helpful here. 
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m 



An algorithm which runs directly in the chiral limit will be very 
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